Simulation apparatus, simulation method, and non-transitory computer readable medium storing program

ABSTRACT

An input unit that receives shape definition data that defines a shape of a simulation target object and a mesh division condition relating to a mesh size when a mesh division of the simulation target object is performed. A processing unit performs the mesh division of the simulation target object based on the shape definition data and the mesh division condition which are input to the input unit and disposes a particle at a node of a generated mesh to perform time evolution of a state of the particle based on an interaction potential between the particles. The processing unit varies a time step width for the time evolution of the state of the particle according to a distance between the particles.

RELATED APPLICATIONS

The content of Japanese Patent Application No. 2018-231799, on the basisof which priority benefits are claimed in an accompanying applicationdata sheet, is in its entirety incorporated herein by reference.

BACKGROUND Technical Field

Certain embodiment of the present invention relates to a simulationapparatus, a simulation method, and a non-transitory computer readablemedium storing a program.

Description of Related Art

In mechanical structure analysis using dynamic explicit methods such asmolecular dynamics method and renormalization group molecular dynamicsmethod, it is preferable to perform time evolution with a time stepwidth or less according to a mesh size (interparticle distance) in ordernot to diverge a calculation. The smaller the mesh size (the shorter theinterparticle distance), the smaller a preferred time step width. Whenthe mesh size of a system to be analyzed differs depending on alocation, the time step width is required to be reduced according to thesmallest mesh size in the system. For the reason, a calculation amountincreases and an analysis time becomes longer.

The related art discloses a living body simulation apparatus thatperforms analysis by coupling molecular motion in a sarcomere and amuscle continuum model. A time step width in the Monte Carlo simulationprocessing of the sarcomere model is different from a time step width ina finite element analysis of the continuum model.

SUMMARY

According to an embodiment of the present invention, there is provided asimulation apparatus including: an input unit that receives shapedefinition data that defines a shape of a simulation target object and amesh division condition relating to a mesh size when a mesh division ofthe simulation target object is performed; and a processing unit thatperforms the mesh division of the simulation target object based on theshape definition data and the mesh division condition which are input tothe input unit and disposes a particle at a node of a generated mesh toperform time evolution of a state of the particle based on aninteraction potential between the particles. The processing unit variesa time step width for the time evolution of the state of the particleaccording to a distance between the particles.

According to another aspect of the invention, there is provided asimulation method including: performing a mesh division of a shape of asimulation target object and disposing a particle at a node of agenerated mesh; and varying a time step width for time evolution of astate of the particle according to a distance between the particles toperform the time evolution of the state of the particle based on aninteraction potential between the particles.

According to yet another aspect of the invention, there is provided anon-transitory computer readable medium storing a program that causes acomputer to execute a process. The process includes a function ofacquiring shape definition data that defines a shape of a simulationtarget object and a mesh division condition relating to a mesh size whena mesh division of the simulation target object is performed; a functionof performing the mesh division of the simulation target object based onthe shape definition data and the mesh division condition which areinput to the input unit and disposing a particle at a node of agenerated mesh to perform time evolution of a state of the particlebased on an interaction potential between the particles; and a functionof varying a time step width for the time evolution of the state of theparticle according to a distance between the particles.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram of a simulation apparatus according to anembodiment.

FIG. 2 is a flowchart of a simulation method according to theembodiment.

FIG. 3 is a graph showing an example of a relationship between a naturalperiod P and a time step width Δt.

FIG. 4 is a timing chart for executing a first loop L1 to a fourth loopL4.

FIG. 5 is a diagram showing an example of a plurality of particlesdisposed at a plurality of nodes of a mesh generated from a shape of asimulation target object.

FIG. 6 is a flowchart showing a detailed procedure of step S6 (FIG. 2).

FIG. 7 is a flowchart showing a procedure of processing of the firstloop L1 (step S67 in FIG. 6).

FIG. 8A is a flowchart showing a processing procedure of step S66 (FIG.6), and FIG. 8B is a flowchart showing a processing procedure of stepS64 (FIG. 6).

FIG. 9 is a flowchart showing a processing procedure of step S62 (FIG.6).

FIG. 10A is a diagram showing shape models of two objects to besimulated, and FIG. 10B is a diagram showing shape models in a statewhere the two objects are in contact with each other.

FIG. 11A is a schematic perspective view of a verification model, andFIG. 11B is a graph showing a simulation result of a temporal change ina position of an unfixed end portion of the verification model.

FIG. 12 is a graph showing calculation times per one time step requiredfor an intra-object interaction calculation and a position and velocitycalculation in comparison between the embodiment and a comparativeexample.

DETAILED DESCRIPTION

In the example of the related art, the time step width in the finiteelement analysis of the continuum model using the dynamic explicitmethod is constant regardless of the mesh size of the continuum model.For the reason, it is required to reduce the time step widthcorresponding to the smallest mesh size, and thus it is difficult toreduce the analysis time.

There is a need for providing a simulation apparatus, a simulationmethod, and a non-transitory computer readable medium storing a programcapable of reducing an analysis time.

Next, a simulation apparatus and a simulation method according to anembodiment will be described with reference to FIGS. 1 to 12.

FIG. 1 is a block diagram of the simulation apparatus according to theembodiment. The simulation apparatus according to the embodimentincludes an input unit 20, a processing unit 21, an output unit 22, anda storage unit 23. A simulation condition and the like are input fromthe input unit 20 to the processing unit 21. Further, various commandsare input from an operator to the input unit 20. The input unit 20 iscomposed of, for example, a communication apparatus, a removable mediareading apparatus, and a keyboard.

The processing unit 21 performs a mesh division of a shape of asimulation target object based on an input simulation condition andcommand and virtually disposes a particle at a node of a mesh to performa simulation by a molecular dynamics method. Further, a simulationresult is output to the output unit 22. The simulation result includesinformation representing a temporal change in the shape of thesimulation target object. The processing unit 21 includes, for example,a computer, and a program for causing the computer to execute thesimulation by the molecular dynamics method is stored in the storageunit 23. The output unit 22 includes a communication apparatus, aremovable media writing apparatus, a display, and the like.

FIG. 2 is a flowchart of the simulation method according to theembodiment. The processing of each step shown in FIG. 2 is realized bythe processing unit 21 (FIG. 1) executing the program stored in thestorage unit 23.

First, the processing unit 21 acquires shape definition data thatdefines the shape of the simulation target object, physical propertyvalues (physical property values representing density, degree ofelasticity, and the like) of the simulation target object, an initialcondition, boundary condition, mesh division condition, and the like ofthe simulation (step S1). These pieces of data are input to the inputunit 20 and are acquired by the processing unit 21 from the input unit20. When the shape definition data is acquired, the processing unit 21performs the mesh division of the shape defined by the shape definitiondata based on a mesh division condition to generate a shape model (stepS2). For example, a known mesh division algorithm can be used for themesh division processing. In the embodiment, the simulation targetobject is divided into tetrahedral meshes.

When the mesh is generated, the processing unit 21 virtually disposesthe particle at the node of the mesh, assigns mass to each particle, anddetermines an interaction potential between the particles (step S3). Themass of the particle is determined such that the density of the shapemodel is the same as the density of the simulation target object basedon, for example, the density of the simulation target object and theparticle distribution. For example, when the density of the simulationtarget object is uniform, the mass of the particle is relatively smallin a region where the particle distribution density is high, and themass of the particle is relatively large in a region where the particledistribution density is low.

For example, the potential of a spring-mass model is used as theinteraction potential between the particles. A method of determining aspring constant is described in detail in, for example, JapaneseUnexamined Patent Publication No. 2009-37334. Here, the method ofdetermining the spring constant will be briefly described.

First, a Voronoi polyhedron analysis is performed on the particledisposed at the node of the tetrahedral mesh. An area of an interfacecrossing a line segment with particles i and j as both ends of aplurality of interfaces constituting the Voronoi polyhedron isrepresented by S_(ij). It is possible to determine a spring constantk_(ij) from a Young's modulus and the area S_(ij) of the simulationtarget object.

Next, a period of a natural vibration (natural period) is calculated foreach of two particles (particle pair) connected to each other by aspring (step S4). In general, since the pieces of mass of two particlesconstituting a particle pair are not the same, the natural period in astate where one particle is fixed is different from the natural periodin a state where the other particle is fixed. For example, a shorternatural period value of the two natural periods is employed as thenatural period of the particle pair.

Based on the natural period calculated in step S4, the processing unit21 determines, for each particle pair, a time step width when timeevolution of a state of the particle is calculated (step S5). Theshorter the natural period, the more the state of the particle changesin a constant time. In order to analyze the change in the state of theparticle with high accuracy, the time step width may be reduced as thenatural period is shorter. For example, the time step width may beselected from a range of 1/20 or more to 1/10 or less of the naturalperiod. A plurality of time step widths are defined for one shape model.

When the time step width is determined, the processing unit 21 analyzesthe change in the state of the particle with the determined time stepwidth (step S6). Specifically, an equation of motion is solved to obtainthe temporal change in the velocity and position of the particles. Whenthe analysis ends, the analysis result is output to the output unit 22(FIG. 1) (step S7).

FIG. 3 is a graph showing an example of a relationship between a naturalperiod P and a time step width Δt. The horizontal axis represents thenatural period P, and the vertical axis represents the time step widthΔt. For example, when the natural period P is less than P1, equal to orlarger than P1 and less than P2, equal to or larger than P2 and lessthan P3, and equal to or larger than P3, the time step width Δt isrespectively set to Δt₁, Δt₂, Δt₃, and Δt₄. Here, a magnituderelationship of Δt₁<Δt₂<Δt₃<Δt₄ is established. Further, Δt₄ is aninteger multiple of Δt₃, Δt₃ is an integer multiple of Δt₂, and Δt₂ isan integer multiple of Δt₁.

For a particle pair having a natural period P of P3 or more, processingof a fourth loop L4 is executed with the time step width Δt₄ to performthe calculation of a force acting between the particles and the timeevolution of the velocity. For a particle pair having the natural periodP of P2 or more and less than P3, processing of a third loop L3 isexecuted with the time step width Δt₃. For a particle pair having thenatural period P of P1 or more and less than P2, processing of a secondloop L2 is executed with the time step width Δt₂. For a particle pairwhose natural period P is less than P1, processing of a first loop L1 isexecuted with the time step width Δt₁.

FIG. 4 is a timing chart for executing the first loop L1 to the fourthloop L4. The processing of the first loop L1 is executed for each timestep width Δt₁. In other words, when the processing of the first loop L1is executed once, the time evolution is performed by the time step widthΔt₁. The processing of the second loop L2 is executed every time theaccumulation of the time evolution increases by the time step width Δt₂.The processing of the third loop L3 is executed every time theaccumulation of the time evolution increases by the time step width Δt₃.The processing of the fourth loop L4 is executed every time theaccumulation of time evolution increases by the time step width Δt₄. Onetime of processing of performing the time evolution by the time stepwidth Δt₁ is referred to as one time step.

FIG. 5 is a diagram showing an example of a plurality of particlesdisposed at a plurality of nodes of the mesh generated from the shape ofthe simulation target object. A plurality of particle pairs are definedfrom the plurality of particles shown in FIG. 5, and any one of thefirst loop L1 to the fourth loop L4 is associated with each of theplurality of particle pairs based on the natural period. For example,the first loop L1 is associated with particle pairs AB and BC. Thesecond loop L2 is associated with particle pairs BD, BE, and CE. Thethird loop L3 is associated with particle pairs DE and DF. The fourthloop L4 is associated with particle pairs EF and EG.

In the time step in which only the first loop L1 (FIG. 4) is executed,the time evolution of the velocity of the particles and the calculationof the force are performed for the particle pairs AB and BC. In the timestep in which the first loop L1 and the second loop L2 (FIG. 4) areexecuted, the time evolution of the velocity of the particles and thecalculation of the force are also performed for the particle pairs BD,BE, and CE. In the time step for executing the first loop L1, the secondloop L2, and the third loop L3 (FIG. 4), the time evolution of thevelocity of the particles and the calculation of the force are alsoperformed for the particle pairs DE and DF. In the time step in whichall the first loop L1 to the fourth loop L4 are executed, the timeevolution of the velocity of the particles and the calculation of theforce are also performed for the particle pairs EF and EG.

FIG. 6 is a flowchart showing a detailed procedure of the processing ofanalyzing the change in the state of the particle (step S6 in FIG. 2).In the following description, FIG. 5 will be referred to as appropriate.

The processing of the fourth loop L4 is executed until a calculation endcondition of the simulation is satisfied. As the calculation endcondition, for example, the number of time steps, the calculation time,and the like are provided. When the calculation end condition is notsatisfied, first, the equation of motion is solved for the particles E,F, and G included in the particle pairs EF and EG (FIG. 5) to beprocessed in the fourth loop L4 to perform the time evolution of thevelocity of each particle by half of the time step width Δt₄ (step S61).

After step S61 is executed, the equation of motion is solved for theparticles D, E, and F included in the particle pairs DE and DF (FIG. 5)to be processed in the third loop L3 to perform the time evolution ofthe velocity of each particle by half of the time step width Δt₃ (stepS63). Here, the time evolution of the latest velocity after the timeevolution is further performed for the particles on which the timeevolution of velocity is performed in step S61, for example, theparticles E and F. Similarly, the time evolution of the latest velocityis further performed also in the subsequent processing of the timeevolution of the velocity.

After step S63 is executed, the equation of motion is solved for theparticles B, C, D, and E included in the particle pairs BD, BE, and CE(FIG. 5) to be processed in the second loop L2 to perform the timeevolution of the velocity of each particle by half of the time stepwidth Δt₂ (step S65).

After step S65 is executed, the processing of the first loop L1 isexecuted for all the particles (step S67). Details of the processing ofthe first loop L1 will be described below with reference to FIG. 7.

When the processing (time step) of the first loop L1 is repeated aplurality of times and an accumulated value of the time evolution afterstep S65 increases by Δt₂, post-processing of the second loop L2 isexecuted (step S66). In the post-processing of the second loop L2 (stepS66), the time evolution of the velocity is further performed by half ofthe time step width Δt₂. The post-processing of the second loop L2 (stepS66) will be described in detail below with reference to FIG. 8A. Asdescribed above, the method of executing the time evolution of thevelocity by half of the time step width is referred to as a velocityVerlet method.

When the processing of the second loop L2 is repeated a plurality oftimes and the accumulated value of the time evolution after step S63increases by Δt₃, the post-processing of the third loop L3 is executed(step S64). In the post-processing of the third loop L3 (step S64), thetime evolution of the velocity is further performed by half of the timestep width Δt₃. The post-processing of the third loop L3 (step S64) willbe described in detail below with reference to FIG. 8B.

After the processing of the third loop L3 is repeated a plurality oftimes, when the accumulated value of the time evolution after step S61increases by Δt₄, the post-processing of the fourth loop L4 is executed(step S62). In the post-processing of the fourth loop L4 (step S62), thetime evolution of the velocity is further performed by half of the timestep width Δt₄. The post-processing of the fourth loop L4 (step S62)will be described in detail below with reference to FIG. 9.

FIG. 7 is a flowchart showing a procedure of the processing of the firstloop L1 (step S67 in FIG. 6). First, the equation of motion is solvedfor the particles A, B, and C included in the particle pairs AB and BC(FIG. 5) to be processed in the first loop L1 to perform the timeevolution of the velocity of each particle by half of the time stepwidth Δt₁ (step S671). Next, the time evolution of the positions of allthe particles is performed by the time step width Δt₁ based on thelatest velocity of each particle (step S672).

After the time evolution of the positions of all the particles isperformed, determination is made whether there is a contact betweenobjects (step S673). For example, when the simulation target object iscomposed of a plurality of objects, the force based on the interactionpotential determined in step S3 (FIG. 2) acts between particles in thesame object. The force does not act between particles of differentobjects in a state where the objects are not in contact with each other.However, when two objects are in contact with each other, the force isexerted between particles positioned at a contact point of the objects.When determination is made in step S673 that there is a contact betweenthe objects, in the next time step, the interaction potential actingbetween particles of different objects is introduced to perform the timeevolution of the velocity and position for the particle positioned atthe contact point of the objects. The determination of contact betweenthe objects and the interaction potential between the particles ofdifferent objects will be described below with reference to FIGS. 11Aand 11B.

Next, the force acting on the particles is calculated based on theposition of the particle after the time evolution (step S674). In thiscase, the force generated based on the interaction potential between theparticles is calculated only for the particle pairs AB and BC (FIG. 5)to be processed in the first loop L1. Further, for a particle on whichexternal forces such as gravity and a load provided as a boundarycondition acts, these external forces are calculated. For the particlepairs to be processed in the second loop L2 to the fourth loop L4, theforce based on the interaction potential is not calculated. When theequation of motion is solved for each particle, the calculation isperformed based on a combined force obtained by combining the forcesacting on the particles.

When the force acting on the particles is obtained, the equation ofmotion is solved for each of the particles A, B, and C constituting theparticle pairs AB and BC to be processed in the first loop L1 to performthe time evolution of the velocity of the particles by half of the timestep width Δt₁ (step S675). In each step of S671 and S675, the timeevolution of the velocity of the particles is performed by half of thetime step width Δt₁. Therefore, the time evolution of the velocity ofthe particles is performed by the time step width Δt₁.

FIG. 8A is a flowchart showing a processing procedure of step S66 (FIG.6). First, the force acting on the particles is calculated (step S661).In this case, the force generated based on the interaction potentialbetween the particles is calculated only for the particle pairs BD, BE,and CE (FIG. 5) to be processed in the second loop L2. When the forceacting on the particles is obtained, the equation of motion is solvedfor each of the particles B, C, D, and E constituting the particle pairsBD, BE, and CE (FIG. 5) to be processed in the second loop L2 to performthe time evolution of the velocity of the particles by half of the timestep width Δt₂ (step S662). In two steps of S65 (FIG. 6) and S662, thetime evolution of the velocity of the particles is performed by the timestep width Δt₂.

FIG. 8B is a flowchart showing a processing procedure of step S64 (FIG.6). Similarly to the case of step S66, the force is calculated for theparticle pairs DE and DF (FIG. 5) to be processed in the third loop L3(step S641). Thereafter, the equation of motion is solved for each ofthe particles D, E, and F constituting the particle pairs DE and DF(FIG. 5) to be processed in the third loop L3 to perform the timeevolution of the velocity of the particles by half of the time stepwidth Δt₃ (step S642). In two steps of S63 (FIG. 6) and S642, the timeevolution of the velocity of the particles is performed by the time stepwidth Δt₃.

FIG. 9 is a flowchart showing a processing procedure of step S62 (FIG.6). Similar to the case of step S66, the force is calculated for theparticle pairs EF and EG (FIG. 5) to be processed in the fourth loop L4(step S621). Thereafter, the equation of motion is solved for each ofthe particles E, F, and G constituting the particle pairs EF and EG(FIG. 5) to be processed in the fourth loop L4 to perform the timeevolution of the velocity of the particles by half of the time stepwidth Δt₄ (step S622). In two steps of S61 (FIG. 6) and S622, the timeevolution of the velocity of the particles is performed by the time stepwidth Δt₄.

After step S622, the energy of a system to be simulated is calculated(step S623). The energy includes kinetic energy of the particles,elastic energy of the particle pairs, contact energy between theobjects, and input and output energy which is input and output betweenthe outside and the system. Due to the law of conservation of energy, asum of the pieces of energy is substantially constant. Conversely, whenthe sum of the pieces of energy is substantially constant, thesimulation result is assumed to be normal.

After step S623, the simulation result is output to the output unit 22(FIG. 1). The information output to the output unit 22 may include, forexample, a temporal change in the shape of the simulation target object,a temporal change in a position of a representative particle, a temporalchange in energy obtained in step S623, and the like.

Next, the method of determining the presence or absence of contactbetween the objects (step S673 in FIG. 7) will be described withreference to FIGS. 10A and 10B.

FIG. 10A is a diagram showing shape models of two objects to besimulated. The shape model of one object 30 is represented by anaggregate of a plurality of particles 31, and the shape model of theother object 40 is represented by an aggregate of a plurality ofparticles 41. When the object 30 and the object 40 are not in contactwith each other, the particles 31 constituting the one object 30 and theparticles 41 constituting the other object 40 do not interact with eachother.

Whether the two objects 30 and 40 are in contact with each other can bedetermined based on, for example, information on a distance between asurface of the one object 30 and a surface of the other object 40.

FIG. 10B is a diagram showing shape models in a state where the object30 and the object 40 are in contact with each other. The distancebetween the surfaces of the two objects 30 and 40 is shorter than in thecase of FIG. 10A.

Next, when the two objects 30 and 40 are in contact with each other, theinteraction between the particles 31 constituting the shape model of theone object 30 and the particles 41 constituting the shape model of theother object 40 will be described.

When the objects 30 and 40 are in contact with each other, a spring-masspotential is adapted as the interaction potential such that a repulsiveforce acts between the particles 31 and 41 near the contact point. Thespring constant is determined by how much a biting amount betweenobjects (calculation error) can be allowed. When the calculation errorcan be allowed, the spring constant can be reduced.

When the time evolution of the velocity of each particle is performed,the interaction potential is considered for the particles 31 and 41 nearthe contact point. The time step width Δt for the time evolution isdetermined based on the natural period in the same manner as theparticle pair in the same object. Determination is made that thecalculation of the time evolution of the velocity and the calculation ofthe force for a pair of the particles 31 and 41 are performed by whichof the first loop L1 to the fourth loop L4 based on the natural period.

Next, an excellent effect of the embodiment will be described. In theabove embodiment, the time step width of time evolution is set accordingto the natural period for each particle pair based on a plurality ofparticles constituting the shape model to be simulated. In a regionwhere the distance between the particles is long, the mass assigned toeach of the particles is large since the distribution density of theparticles is small. For the reason, the natural period becomes long.That is, the natural period of the particle pair depends on the distancebetween the particles. Therefore, it can be said that the setting of thetime step width according to the natural period in the above embodimentis to set the time step width according to the distance between theparticles.

As described above, since the time step width is increased for theparticle pair having a relatively long distance between particles in oneshape model, it is possible to reduce the calculation load as a whole.Further, since the time step width is determined based on the naturalperiod of the particle pair, it is possible to easily determine apreferable time step width as compared with the case where the time stepwidth is determined by trial and error.

In the above embodiment, considering the contact of two objects, thetime step width of the time evolution adapted to the contact point isdetermined by the same method as the time step width adapted to theparticle pair in the object. For this reason, it is possible to reducethe calculation load also when the above embodiment is adapted to asimulation relating to a many-body contact of a mechanism structuresystem.

Next, a simulation performed to verify the effects of the aboveembodiment will be described with reference to FIG. 11A to FIG. 12.

FIG. 11A is a schematic perspective view of a verification model. Atemporal change in the shape of a plate-like member is simulated in astate where one end portion 50 in the longitudinal direction of theplate-like member that is long in one direction is fixed and a load 52in the thickness direction is applied to the other end portion 51. Thatis, a model in which the load 52 is applied to the tip of a cantileverstructure is employed as the verification model. The mesh size isrelatively reduced near the fixed end portion 50. The number of meshnodes (particles) is 15,283, and the number of particle pairs is 93,740.

Three different time step widths are set to perform the time evolutionof the state of the particle. The time step widths are denoted as Δt₁,Δt₂, and Δt₃ in ascending order. The time step widths Δt₂ and Δt₃ arerespectively set to 2 times and 4 times Δt₁. It is possible to adapt thetime step width Δt₃ to about 60% of the particle pairs and the time stepwidth Δt₂ to about 25%. The time step width Δt₁ is adapted to theremaining about 15% of the particle pairs. For comparison, a simulationby a comparative example in which the time step width Δt₁ is adapted toall the particle pairs is also performed.

FIG. 11B is a graph showing a simulation result of the temporal changein a position of the unfixed end portion 51 (FIG. 11A) of theverification model. The horizontal axis represents time in the unit“ms”, and the vertical axis represents the position in the unit “mm”. Inthe graph of FIG. 11B, solid circle symbols and white circle symbolsrespectively indicate the simulation results of the embodiment and thecomparative example. It is confirmed that both of the circle symbolssubstantially overlap each other and the same accuracy as in thecomparative example is obtained in the embodiment.

FIG. 12 is a graph showing the calculation time per one time steprequired for an intra-object interaction calculation and a position andvelocity calculation in comparison between the embodiment and thecomparative example. The intra-object interaction calculationcorresponds to the force calculation performed for each particle pair(step S674 in FIG. 7, step S661 in FIG. 8A, step S641 in FIG. 8B, andstep S621 in FIG. 9). The position and velocity calculation correspondsto the calculations of the time evolution of the velocity and the timeevolution of the position.

In the case of the model shown in FIG. 5 as an example, the particle Eis included in the particle pairs CE and BE corresponding to the secondloop L2, the particle pair DE corresponding to the third loop L3, andthe particle pairs EF and EG corresponding to the fourth loop L4. Forthis reason, the time evolution of the velocity of the particle E isperformed in all the second loop L2, the third loop L3, and the fourthloop L4. Therefore, in the processing of performing the time evolutionof the position and velocity, a calculation time of the embodiment islonger than a calculation time of the comparative example.

In the intra-object interaction calculation, since the time step widthis increased for many particle pairs, the calculation time of theembodiment is reduced to about 40% of the calculation time of thecomparative example. Regarding the entire calculation time including theupdate of the position and velocity, the calculation time of theembodiment is about 46% of the calculation time of the comparativeexample. It is confirmed that the calculation time can be reduced ascompared with the comparative example by adapting the simulation methodaccording to the embodiment.

Next, a modification example of the above embodiment will be described.In the above embodiment, the four kinds of time step widths Δt₁, Δt₂,Δt₃, and Δt₄ are used as the time step width Δt based on the naturalperiod. However, the time step widths may be two kinds or more.

The above embodiment is an example, and the invention is not limited tothe above embodiment. For example, it is apparent to those skilled inthe art that various changes, improvements, combinations, and the likecan be made.

It should be understood that the invention is not limited to theabove-described embodiment, but may be modified into various forms onthe basis of the spirit of the invention. Additionally, themodifications are included in the scope of the invention.

What is claimed is:
 1. A simulation apparatus comprising: an input unitthat receives shape definition data that defines a shape of a simulationtarget object and a mesh division condition relating to a mesh size whena mesh division of the simulation target object is performed; and aprocessing unit that performs the mesh division of the simulation targetobject based on the shape definition data and the mesh divisioncondition which are input to the input unit and disposes a particle at anode of a generated mesh to perform time evolution of a state of theparticle based on an interaction potential between the particles,wherein the processing unit varies a time step width for the timeevolution of the state of the particle according to a distance betweenthe particles.
 2. The simulation apparatus according to claim 1, whereinthe processing unit determines an interaction potential between theparticles and mass of the particle based on a physical property valuerepresenting a degree of elasticity of the simulation target object, adensity of the simulation target object, and a shape of the generatedmesh, and determines the time step width for each particle pair based ona period of a natural vibration period obtained from the interactionpotential between the particles and the mass of the particle.
 3. Thesimulation apparatus according to claim 1, wherein the simulation targetobject includes a plurality of objects, and wherein the processing unitdetermines whether there is a contact between the objects by timeevolution, and when determination is made that there is a contactbetween the objects, determines the time step width for the particlepositioned at a contact point based on a period of a natural vibrationobtained from mass of the particle disposed at the contact point and aninteraction potential between the particles.
 4. A simulation methodcomprising: performing a mesh division of a shape of a simulation targetobject and disposing a particle at a node of a generated mesh; andvarying a time step width for time evolution of a state of the particleaccording to a distance between the particles to perform the timeevolution of the state of the particle based on an interaction potentialbetween the particles.
 5. A non-transitory computer readable mediumstoring a program that causes a computer to execute a function ofacquiring shape definition data that defines a shape of a simulationtarget object and a mesh division condition relating to a mesh size whena mesh division of the simulation target object is performed; a functionof performing the mesh division of the simulation target object based onthe shape definition data and the mesh division condition which areinput to the input unit and disposing a particle at a node of agenerated mesh to perform time evolution of a state of the particlebased on an interaction potential between the particles; and a functionof varying a time step width for the time evolution of the state of theparticle according to a distance between the particles.